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Current-free double layers of the type reported in plasmas in the presence of an expanding magnetic field [C. 
Charles and R. W. Boswell, Appl. Phys. Lett. 82, 1356 (2003)] are modeled theoretically and with particle- 
in-cell/Monte Carlo simulations. Emphasis is placed on determining what mechanisms affect the electron 
velocity distribution function (EVDF) and how the EVDF influences the double layer. A theoretical model 
is developed based on depletion of electrons in certain velocity intervals due to wall losses and repletion of 
these intervals due to ionization and elastic electron scattering. This model is used to predict the range of 
neutral pressures over which a double layer can form and the electrostatic potential drop of the double layer. 
These predictions are shown to compare well with simulation results. 



PACS numbers: 52.27.-h,52.65.Rr,52.75.Di 
I. INTRODUCTION 

Double layers are adjacent regions of net positive 
and negative charge that form distant from the physi- 
cal boundaries of a plasma. They typically provide an 
electrostatic boundary that separates plasmas with dif- 
ferent properties. There are several varieties of double 
layerspEl some of which have been studied since the ear- 
liest days of plasma physics researchP One categoriza- 
tion is that double layers can be either current-carrying 
or current-free. The current-free variety was predicted 
theoretically in the early 1980sp and these were later 
observed experimentallyPS Recently, a renewed interest 
in current-free double layers^^ has arisen in part be- 
cause of their application to electrostatic thrusters for 
spacecraft propulsiorPSHUl and auroral physics P^l 

These recent current-free double layer experiments^^ 
consist of an insulated source chamber connected to a 
larger volume expansion chamber that is metallic and 
grounded; see Fig. [T] An approximately constant axial 
magnetic field is applied to the source chamber, which 
diverges near the boundary between the source and ex- 
pansion chambers. Plasma is generated in the source 
chamber by applying rf waves with an antenna. Current- 
free double layers have been measured in the reg ion of 
divergent magnetic field in this configuration!-^ It has 
also been confirmed that these double layers generate an 
ion beam in the expansion chamber that has a flow s peed 
typically a few times faster than the ion sound speedP^El 

Analytic models of current-free double layers in ex- 
panding plasmas have been proposed by ChenJ^SI Lieber- 
man et a/j^3 Goswami et aZp3 and Ahedo and SanchezP^l 
These are fundamentally different in that each makes a 
different assumption for the electron velocity distribu- 
tion function (EVDF). Chen considers just the upstream 
region and assumes that electrons are Maxwellian.^ 
Lieberman et al consider two populations of electrons 
upstream: a thermal (Maxwellian) population and an 
additional half-Maxwellian beam population.^ The up- 
stream electrons in Goswami et al are counter-streaming 



Maxwellian beamsP^ Ahedo and Sanchez assume a two- 
temperature Maxwellian distribution characterized by 
hot and cold populations!^ Double layer formation is 
sensitive to the EVDF, so each of these theories predicts 
different double layer parameters such as the potential 
drop and resultant ion beam properties. 

An accurate model of the EVDF, and experimental 
verification of it, is needed to provide a foundation for 
a comprehensive analytic model of the experiments. In 
particular, verification of the electron beams assumed to 
be present in the source chamber in Refs. [55] or [57] is 
lacking. Unfortunately, diagnosing the EVDF is diffi- 
cult to do experimentally. Essentially the only diagnos- 
tic available is a Langmuir probe, but this is typically 
limited to measuring the electron energy distribution 
function (EEDF) rather than the EVDF. Another lim- 
itation of Langmuir probes is that current- volt age char- 
acteristics get noisy for energies greater than a couple 
of electron temperatures. Previous measurements have 
given ambiguous results concerning electron beams in the 
source region. Early work with Langmuir probes pro- 
vided some "preliminary" evidence of an electron beam 
very close to the sheath of the source chamber!^! Other 
indirect measurements associated with an ionization in- 
stability also appeared to suggest that electron beams 
were present P3 However, more recent Langmuir probe 
measurements have found no evidence of beamsPSMl In- 
stead, these found a Maxwellian EEDF that was depleted 
in density beyond the double layer potential energy. 

In this work, we develop a model for the EVDF in 
an expanding plasma and compare the results with PIC 
simulations. We concentrate on a simplified geometry 
that has only one spatial dimension, but is 3D in veloc- 
ity phase-space. The boundary conditions on the geomet- 
ric dimension are an insulating wall at one end (source 
chamber) and a conducting wall at the other (expansion 
chamber). The analytic model accounts for depletion of 
velocity phase-space intervals due to loss of electrons to 
the boundaries, as well as partial repletion of these inter- 
vals due to ionization sources and scattering. The model 
predicts that a current-free double layer can only exist 
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FIG. 2. Sketch of a typical potential profile for a current-free 
double layer in the experimental geometry shown in Fig. [T] 



FIG. 1. Schematic diagram of a common experimental appa- 
ratus for studying current-free double layers in an expanding 
magnetic field. 



over a finite range of neutral pressures. It can also be 
used to predict the double layer and sheath potential 
drops based on the electron temperature. 



The PIC code, named phoenix, uses the same ID in 
space, 3D in velocity phase-space geometry that the an- 
alytic model is based on. Collision processes are mod- 
eled using a Monte Carlo algorithm and energy is input 
with a method that simulates inductive electron heating 
in a velocity-space direction perpendicular to the geo- 
metric domain. Plasma expansion is modeled by invok- 
ing a loss profile in the downstream region. Aside from 
details of the loss profile used, PHOENIX has been de- 
signed to be identical to the code developed by Meige 
et aZPESI We find that the EVDFs calculated using the 
PIC code differ substantially from those assumed in pre- 
vious literaturel^SHIH Electron beams are not observed, 
which is consistent with the most recent Langmuir probe 
measurements EDEH The EVDFs in the simulations are 
shown to agree with the analytic model based on de- 
pletion due to wall losses and partial repletion due to 
scattering. 



This paper is organized as follows. Section |H| develops 
a model for the EVDF starting from academic examples 
which highlight the physics behind the maximum dou- 
ble layer potential drop, as well as the minimum and 
maximum neutral pressures that can support a double 
layer. A model of the whole ID simulation domain is 
given in Sec. |IID| After describing the phoenix code in 
Sec. |III| the simulation results are provided in Sec. |IV| 
This section also contains a discussion of how the simu- 
lated EVDFs relate to previous work and compare with 
the analytic model of Sec. [II} The results are summarized 
in Sec. N\ 



II. MODEL OF DOUBLE LAYER FORMATION AND 
THE EVDF 

The electrostatic potential profile along the axis of 
previous current-free double layer experiments is shown 
schematically in Fig. [2pISl The boundary on the upstream 
side (region 2) of the experiments is insulating, while the 
downstream boundary of the larger expansion chamber 
is conducting and grounded, which we take to be the ref- 
erence potential. Of course, the experiments, which are 
cylindrical, have radial profiles in the transverse direction 
that can affect the details of the axial profile at differ- 
ent radial positionsF^M] Although consideration of these 
3D affects are necessary in order to quantitatively model 
the experiments, we consider a simplified ID model here. 
Our goal is to identify what mechanisms influence the 
EVDF and, as a result, the double-layer and sheath po- 
tentials. The simulations presented in Sec. IV also use 



the ID geometry and thus provide a proving ground for 
comparison to this model. 

Since the upstream boundary is insulating, it must 
collect equal fluxes of electrons and ions (assumed to 
be singly charged here) during steady-state operation. 
The only other physical boundary in this system is the 
grounded downstream wall, thus it too will collect equal 
electron and ion currents (we assume no external elec- 
tron or ion sources, the only source is ionization which 
produces electrons and ions in equal numbers). In the 
absence of any current sources or sinks, a consequence of 
the current-free boundary conditions is that the double 
layer must also be current-free. If the EVDF [f e ,x(v x )] 
is known at the positions P = S2, DL2 and SI, denot- 
ing the upstream sheath edge, the upstream double-layer 
edge, and the downstream sheath edge respectively, the 
current-free condition 



dv x 



fe,x (•£ — P, 



= e 1/2 n P c s , P 



(1) 



can be used at each of these locations to determine the 
upstream sheath potential drop (A0 S 2), the double layer 
potential drop (A^dl) an d the downstream sheath po- 
tential drop (Acf> s i). The right side of Eq. is the Bohm 
flux for ions and c s = yjT e /Mi the ion sound speed. 
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FIG. 3. Sketch of the natura l log of the EVDF for the semi- 
infinite domains of Sec. |II A] at positions (a) just upstream of 
the double layer (P=DL2), and (b) just downstream of the 
double layer (P=DL1). 



The e -1//2 term is from the density drop caused by the 
presheath. 

In the following four sections, we develop a model for 
f e<x that can be used in Eq. (fTl) to calculate the dou- 
ble layer potential drop as well as determine a neutral 
pressure range that can support it. Section |II A| starts 
with a simplified geometry in which both regions 1 and 
2 are semi-infinite domains. This geometry, which has 
also been studied by ChenpSl provides a maximum dou- 
ble layer potential drop. The following sections, |II B| and 
|II C[ account for the upstream and downstream walls, 
which leads to predictions for the minimum neutral pres- 
sure and maximum neutral pressure that can support a 
double layer. Section |IID| puts these geometries together 
to form a comprehensive model for the EVDF that ac- 
counts for both upstream and downstream boundaries. 



rected energy to traverse the double layer and escape 
downstream. Downstream, these electrons create a half- 
Maxwellian distribution. The EVDF at position DL2 can 
thus be written 



e,DL2 — 



/TTV Te 



0, 

?VDL2, 



Vx < -VDL 



(2) 



Here ro aj DL2 is a density variable corresponding to the 
— ^dl < v x < oo region of velocity space. If the distri- 
bution were Maxwellian for all velocities, n a DL2 would 
equal the total density t£dl2 = dv x f e ,DL2- However, 
since / e , D L2 = for v x < -u DL) n DL2 < n a ,DL2- Like- 
wise, T e is not equal to the total temperature defined 
from a velocity-space moment of f e , but the two are ap- 
proximately the same as long as e|A0DiJ ^ T e . 

Putting Eq. ^ into the current-free condition of 
Eq. ([I]) provides an expression relating the double layer 
potential drop and the electron temperature 



\n a v L2 v e e-^^ =n DL2 e- 1 / 2 J?±. (3) 



Here v e = \/8T e /(7rm e ) is an average electron speed. 
Solving Eq. M for |A0 DL | yields 



\^\ = --H—c- 1/2 J 2 -^)- w 



n a,DL2 



Recall that 7i a ,DL2 > ?iDL2, but from the definition 

"-DL2 = fToo dv x / e ,DL2 : 



A. Semi-infinite domains: |A0 D l max 

We start with perhaps the simplest conceptual current- 
free double layer configuration. Here it is assumed 
that plasma is generated in an upstream source region 
that is sufficiently large that the plasma has a nom- 
inal Maxwellian distribution (i.e., the source chamber 
is longer than either the electron-electron or electron- 
neutral collision length). The downstream region is as- 
sumed to be infinite and collisionless, so all particles that 
escape the source region remain downstream. This con- 
figuration may be relevant to a thruster operating in 
space, where the thruster is the source and downstream 
is space vacuum. 

The expected EVDF just upstream and downstream of 
the double layer is shown in Fig. [3] for this configuration. 
Here £ x = \m e v x \v x \ is an energy variable that accounts 
for the particle direction. At position DL2 (just upstream 
of the double layer), the distribution is Maxwellian in 
the velocity interval v x > 0, which consists of thermal 
electrons migrating from the upstream region. It is also 
Maxwellian in the interval — i>dl = — v/2e] A^DLl/we < 
v x < 0, which consists of thermal electrons from the 
source that were subsequently reflected from the dou- 
ble layer electric field. The distribution is empty in the 
interval v x < — i>dl since these electrons had enough di- 



n DL2 
?VDL2 



-erfcA 



'e|A</> 



DL 



(5) 



since erfc( v / e| A0 DL |/T e ) - 0{^/m e /Mi) < 1. Thus, 
the double layer potential drop is approximately the 
floating potential of a planar probe 



|A0 



DL 



Te 

2e 



1 + ln 



(6) 



Equation ([6j has previously been derived by CherPSl 
in the context of current-free double layers. Although 
the semi-infinite domain approximation may be useful 
for a thruster operating in space, it is unable to capture 
some features of finite laboratory experiments. Equa- 
tion ([6| provides a maximum potential drop that might 
be expected in the laboratory. Accounting for plasma 
in a finite downstream expansion chamber leads to some 
electrons migrating up the double layer and being accel- 
erated into the source region. These electrons fill in part 
of the otherwise truncated tail of the EVDF. To preserve 
the current balance in this situation, the double layer po- 
tential must be reduced in comparison to Eq. ([6| so extra 
electrons are allowed to leak downstream to balance those 
coming upstream. This effect will be discussed in more 
detail in Sec. iHCl 
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FIG. 4. Sketch of the natural log of the EVDF for the fi- 
nite source region of Sec. II B at (a) the sheath edge of the 
source wall (P=s2), and (b) just upstream of the double layer 
(P=DL2). 



FIG. 5. Sketch of the nat ural l og of the EVDF for the finite 
expansion chamber of Sec. |II C] at positions (a) just upstream 
of the double layer (P=DL2), and (b) just downstream of the 
double layer (P=DL1). 



B. Upstream wall effects: p n 



Next, we consider a geometry with the same semi- 



infinite and collisionless downstream region as Sec. II A 



but allow for a source chamber of finite length. For this 
case, we model the EVDF at position P as 



fx,P = 



IK* 



/TTV Te 




V x < -v DL 
-u DL <v x < v s2 

V S 2 < V x 



(7) 



in which i> s2 = -\/2e|A0 S 2|/m e . The distribution of 
Eq. ([7]) is shown schematically for positions P=s2 and 
P=DL2 in Fig. |] The distribut ion is Maxwellian with 
density n a in the velocity-space interval where particles 
are confined: — vdl < v x < v s2- Outside of this inter- 
val (in the tails) the EVDF is depleted from the nom- 
inal Maxwellian distribution due to losses to the wall 
through the upstream sheath, or to the downstream vac- 
uum through the double layer. These regions get repleted 
in the source primarily due to elastic collisions from the 
perpendicular to parallel direction. Equation Q models 
these tail regions by assigning a different density (rib or 
n c ) to the tail regions. It is assumed that these regions 
can be described by the same temperature as the bulk 
interval. We also assume that the upstream sheath and 
double layer are sufficiently thin that they are approxi- 
mately collisionless. Thus, n CjS 2 = 1.&.DL2 ~ 0. 

With the assumed boundary conditions, the source 
chamber is essentially a plane symmetric discharge. Due 
to this symmetry = w c ,DL2j which implies A^dl = 
Acf> S 2- Applying these assumptions, and putting Eq. |7]) 
into Eq. Q, yields 



T, 



|A0 DL | = - — In 

e \n-cDL2 



rt DL2 _i/2 



27rm c 



(8) 



Aside from the density ratio, «.DL2/ rt c.DL2: Eq. |8]) is sim- 
ply the floating potential from Eq. Q. However, since 
"DL2/n c ,DL2 ~ «-a,DL2/«-c,DL2 > 1, the extra term acts 
to reduce the double layer potential. Equation ^ has a 
viable solution only if 



< :^L e -V2 

™c,DL2 



27TTO e 



< 1. 



(9) 



Equation ([£]) shows that if the there is not enough scat- 
tering in the source region, the discharge cannot be 
maintained. Scattering in the source causes the oth- 
erwise missing tails of the EVDF to be filled in, so 
n c ,DL2/^DL2 = f{K-n/L s ) in which A e _„ is the electron- 
neutral scattering length and L s is the length of the 
source region. The particular functional dependence of 
this relationship depends on details of the scattering cross 
sections. However, if we assume that it has a simple lin- 
ear dependence 



™c,DL2 
n DL2 



L s /X e 
1, 



(10) 



this can be used to estimate the minimum neutral pres- 
sure required to maintain the discharge. Using A e _„ = 
l/(n„<7 e _„) and n n — n a p, in which n = 3.3 x 10 19 [m -3 
mTorr -1 ] and p is in mTorr, Eq. ([9| implies 



Pn 



flo&e—n-Ls 



For neutral pressures less than Eq. (11 
double layer is not predicted to be a steady- 



(11) 

, a current free 
state solution. 



C. Downstream wall effects: p m ax 

If the expansion chamber downstream is finite in ex- 
tent, the sheath at the downstream wall will reflect a 
population of electrons that can migrate back to the dou- 
ble layer. These are subsequently accelerated into the 
source chamber. In addition, scattering in the down- 
stream region can partially replete the velocity space 
interval beyond the downstream sheath cut-off: v s i = 
y/2e\A(j) s i\/m e . The EVDF just up and downstream of 
the double layer is shown in Fig. [5] for this case. At the 
upstream position, the EVDF takes the form 



x,DL2 



y/lTVTe 



«d,DL2- 
n a ,DL2: 



V x < -«DL+sl 
-VDL+sl < V x 



, (12) 



in which vuh+si = \/2e(\A<j) DL \ + |A0 s i|)/m e . The 
EVDF just downstream has the same form, but with 
«DL+si replaced by u al . 
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Putting the EVDF from Eq. ( 12 ) into the current-free 
condition of Eq. yields 



T e l J ny lj2 e- 1 l 2 ^2imi e /M i 

6 V "a,DL2 - n-d,DL2 



(13) 



in which \(f)2\ = |A</>dl| + |A0 s i|. Equation (13) has a 
viable solution only if 



< 



™DL2 



n a,DL2 — 71(i,DL2 



Assuming n ajD L2 ~ ^dl2, Eq. (14) requires 



™<1DL2 
™DL2 



< 1 



-1/2 



2nm e 
Mi 



1. 



(14) 



(15) 



As in Sec. |IIB[ the precise functional dependence of 
^<i,DL2/^DL2 due to scattering in the downstream region 
is difficult to determine. We again assume a simple linear 
form 



"-DL2 



, f Ld/Xe-n, La < A, 

' \ 1, L d > A, 



i e,n 



(16) 



in which is the length of the downstream region. Ap- 
plying the relations A e _ n = l/(n„er e _ n ) and n n — n Q p, 
in which n Q = 3.3 x 10 19 [m~ 3 mTorr -1 ] and p is in 



mTorr, Eqs. (15) and (16) provide an estimate for the 



maximum neutral pressure the current-free double layer 
solution can support 



Pn 



1 

Wnd p — r}.Lfi 



(17) 



Equation (17) shows that when too many electrons mi- 



grate up the double layer from downstream, the double 
layer potential cannot adjust enough to preserve current 
balance. The physics justification for Eqs. (11) and (17) 
are the same as those determining the p min and p max in 
Ref. However, the analysis is different since Ref. US] 
is based on a 3D fluid model which is diffusion domi- 
nated, while this is a ID kinetic model where collisions 



are modeled with the simple linear estimates of Eqs. ( 10 ) 



or (161 



D. Finite ID domain 

The full simulation domain has two boundaries and 
the length of both the source and downstream domains 
can be comparable to A e _„ (depending on the neutral 
pressure). For low neutral pressures, we expect that the 
EVDF will reflect features of losses to both walls in the 
manner depicted in Fig. [6j Figure [6] shows a sketch of 
the expected EVDF at four locations in the simulation 
domain: the upstream sheath edge (s2), just upstream of 
the double layer (DL2), just downstream of the double 
layer (DL1), and the downstream sheath edge (si). As 
the figure demonstrates, this model of depletion due to 
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FIG. 6. Sketch of the n atural log of the EVDF for the finite 
ID domain of Sec. II D at (a) the sheath edge of the source 
wall (P=s2), (b) just upstream of the double layer (P=DL2), 
(c) just downstream of the double layer (P=DL1), and (c) the 
sheath edge of the expansion chamber wall (P=sl). 



wall losses and repletion due to scattering predicts sev- 
eral features of the EVDF that can be tested in the stim- 
ulations. At low neutral pressures, particularly, these 
features should be clearly visible and their location in 
velocity-space can be compared with the predicted val- 
ues dependent on the sheath and double layer potential 
drops. As the neutral pressure is increased, repletion 
becomes more prevalent and velocity-space intervals af- 
fected by wall losses are more quickly filled in. At higher 
neutral pressures, it is expected that the depleted inter- 
vals become more difficult to distinguish until finally the 
downstream region becomes too collisional to support the 
current-free double layer solution. 

The sheath and double layer potential drops can be 
written in terms of the densities of the various intervals 
in velocity space, in a similar manner to Sees. |II B| and 
|IIC[ but the extra velocity-space intervals significantly 
complicate the analysis. The only qualitative difference 
to the analysis of the previous two sections is that ac- 
counting for migration of a small current of downstream 
electrons into the upstream region leads to a slight asym- 
metry in the source region (so |A0 S 2| ~ |A0dl2|, instead 
of |A0 s2 | = \A<f> DL2 \)- We expect that Eqs. (Ill]) and (fl7} 



remain good approximations for the minimum and maxi- 
mum pressure limits, and that the double layer potential 
drop remains close to the floating potential of Eq. ^ 
for intermediate pressures. These estimates will be com- 
pared with simulation data in Sec. |IV| 



III. DESCRIPTION OF THE PHOENIX CODE 



phoenix is a PIC-MCC code that is ID in space and 



G 



Quantity 



Value 



Neutral pressure 

Domain length 

Number of grid cells 

Time step 

Total run time 

Antenna frequency (uj /2tt) 

Antenna current density amplitude 

*7factor 



1 mTorr 
10 cm 
250 

5 x lO^ 11 s 
25 /is 
10 MHz 
100 A/m 2 
8 x 10 8 
1 x 10 6 s~ 



TABLE I. Parameters used for all simulations except those 
shown in Figs. |12| |16| and |17| in which the neutral pressure 
was changed. The gf ac tor was also adjusted for these so that 
the total number of particles in steady-state exceeded 10 . 



3D in velocity phase-space (1D-3V). It is designed to be 
identical to the JanuS code described in Meige et alW^ 
The left wall (source chamber) is a floating boundary, 
which is achieved computationally by inserting a capac- 
itor there. The right wall (expansion chamber) is con- 
ducting, which is implemented by removing all particles 
that reach the cell defining that boundary. Collisions be- 
tween macroparticles (typically representing ~ 10 9 real 
particles) arc simulated with a Monte Carlo technique in- 
cluding the null collision method based on the algorithm 
developed by Vahedi and Surendra.^Sl The gas species 
here is argon. The cross section for electron impact ion- 
ization was taken from Krishnakumar and Sri vast avag™ 
and electron excitation collisions from de Heer et alF^ 
The electron-argon elastic scattering cross sections were 
taken from Ferch et aP^l for 0-20 eV and from de Heer et 
aPZIfor 20-3000 eV. These are also collected in HayashifHI 
The cross sections for argon ion charge-exchange and ion- 
ization collisions are from Phelps.^ 

The plasma is generated by first loading a small num- 
ber of macroparticles (typically 1000) with a spatially 
uniform Maxwellian distribution of temperature 1 eV 
throughout the simulation domain. Electrons are heated 
in a single Cartesian velocity-space direction (y) per- 
pendicular to the spatial dimension (x) using the in- 
ductive heating method described in Meige et aP^ The 
macroparticle density initially increases due to electron- 
neutral ionization collisions. Eventually, a steady-state is 
reached where particle generation balances particle loss. 
This typically occurs within 25 [is and the typical time 
step used is 50 ps. The number of macroparticles in 
steady-state is > 10 5 . The parameters used in all simula- 
tions are summarized in Table ||J except that the neutral 
pressure was varied for the simulations shown in Figs.|12( 



16]and[T7 The factor was also adjusted for these to meet 
the > 10 b macroparticle condition. These calculations 
were performed on a desktop PC, and each run took 2-5 
days. 

In the experiments, a double layer forms due to the 
expansion of the plasma volume. As the volume ex- 
pands, the plasma density drops. If this density drop is 
steep enough, a double layer will form. Since the simula- 



tions have only one spatial dimension, volume expansion 
cannot be simulated self-consistently. Instead, a density 
drop is imposed by removing particles from the system 
at a set frequency defined by a profile and amplitude. In 
Ref. [351 various linear loss profiles were used to generate 
a double layer, but these did not necessarily represent the 
effective loss profile associated with an expanding mag- 
netic field. Here we modify the loss profile to more closely 
resemble a diverging solenoidal magnetic field. 

The vacuum magnetic field on axis from the coil closest 
to the expansion chamber is B [l + (x — x c ) 2 /i? 2 ]~ 3//2 , 
in which R is the coil radius, x c is the axial position of 
the coil and B Q = fi I/(2R) where I is the coil current. 
We assume that the magnetic field is constant inside the 
source chamber, so the magnetic field on axis throughout 
the domain is 

B W= B °{(1 + Xh-^IXTL (18) 

in which X = (x — x c )/R. The volume expansion obeys 
V/V — (r/r ) 2 = B / B fSi so the change in volume sat- 
isfies V^dV/dx = B~ 1 \dB/dx\. Thus, an appropriate 
loss profile for magnetic field expansion has the form 
^loss ~ [v/B a )\dB/dx\, in which v is some characteris- 



tic velocity. For the field of Eq. (18 1, the loss profile is 



Via 



( \ — j 0' ® — x — Xc (-\a\ 

,{X) = { 3v X(l + X 2 )- 5 / 2 , x c < x < L [W) 

in which v Q = v/R. Note that f max ~ 0.86^ o . For all our 
simulations we chose L = 10 cm and x c = L s = 5 cm. In 
the experiments, R/L s ~ 0.17, and in order to preserve 
this ratio we take R = 1.7 cm in the simulations. We 
will also choose v Q = 1 X 10 6 s _1 , which corresponds to 1 
eV electrons (the initial electron temperature). The loss 
profile of Eq. ( 19 ) is shown in Fig. [7J along with the linear 
loss profile used in Ref. 35. Unless otherwise specified, 
the simulation results presented in the following sections 
used the loss profile from Eq. (19 1. 



Although this simplified simulation geometry can pro- 
vide insight into the mechanisms of double layer forma- 
tion and the role of the EVDF, especially in testing the 
model of Sec.|Hj it is not a quantitatively accurate model 
of the experiments. Since the code is ID, it does not cap- 
ture radial effects that have been the topic of recent ex- 
perimental workF^M] Also, the 10 cm length of the sim- 
ulation domain is nearly an order of ma gnit ude shorter 
than the axial length of the experiments!^ Aside from 
these geometrical effects, one also needs to be cognizant 
of the physics limitations of this model when interpreting 
the simulation data. The loss profile is a mock-up of the 
density drop due to an expanding field, but there is no 
actual magnetic field in the simulations. For instance, 
V-B drifts may play a role in the expansion region, but 
are not captured in the simulations. Since the loss pro- 
file removes particles randomly (independent of energy) , 
slower particles are more likely to be removed in the loss 
region. Also, the neutral density is assumed to be uni- 
form and constant, so effects of neutral depletion, which 
may be important in experiments,^ are not captured. 
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FIG. 7. Loss profiles implemented in phoenix to simulate 
plasma volume expansion downstream. The triangular loss 
profile (dashed red line) was used in Ref. (33 and the curve 
representing an expanding magnetic field (solid black line) is 
from Eq. Jl9|. 



IV. SIMULATION RESULTS 

The electrostatic potential and density are shown in 
Fig. [8] for both the linear and expanding magnetic field 
loss profiles from Fig. [JJ The data shown throughout 
this work was averaged over a few rf periods. The den- 
sity and potential profiles are qualitatively similar for 
either loss profile. However, the upstream potential is 
a few volts less for the magnetic field expansion profile 
[from Eq. (19)]. Also, the double layer potential drop 



is steeper and the downstream region more uniform for 
Eq. ( 19 ). These are due to the relative narrowness of the 
magnetic field expansion profile, which is shown in Fig. [7] 
The characteristic step potential profile of a double layer 
is seen in Fig. [8] Figure [9] shows a profile of the charge 
density: p — e(rii — n e ). It has been suggested in pre- 
vious literature that the potential profile of expanding 
plasmas, which are usually deemed "double layers," are 
actually single layers similar to sheathsPS Figure [9] shows 
explicitly adjacent regions of positive and negative space 
charge, which is typically the property used to define a 
double layerP Thus, we conclude that double layers, not 
single layers, are found in these simulations. 



Ion beams and the IVDF 




FIG. 8. Plasma potential and electron density throughout 
the simulation domain. The red dashed line corresponds to 
a simulation that used the linear expansion profile and the 
solid black line to one that used the expanding magnetic field 
profile; see Fig. [7] 




The ion velocity distribution function (IVDF) in the x 
direction is shown as a color map in Fig. [10] through- 
out the simulation domain. In the central source re- 
gion, it is a stationary Maxwellian. Ions are acceler- 
ated by the sheath electric fields at each boundary, so 
the IVDF has a flow shift there and a lower energy tail 
due to ion scattering. A supersonic ion beam is gener- 
ated by the double layer potential drop and this beam 
is maintained at a constant speed downstream (until the 
downstream sheath is reached). The beam speed is ap- 



X [cm] 



FIG. 9. Charge density as a function of position showing the 
adjacent regions of positive and negative space charge that 
define a double layer. The red dashed line corresponds to a 
simulation that used the linear loss profile and the black solid 
line to one that used Eq. ( 19 1. 




x [cm] 



FIG. 10. A color map showing the natural logarithm of the 
ion velocity distribution in the x direction throughout the 
simulation domain. Colors corresponding to higher numbers 
on the color bar represent higher concentration of particles. 
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»« [k.»/s] 

FIG. 11. The ion velocity distribution function in the x direc- 
tion divided by the ion density at the axial locations x =2.5 
cm, 5 cm, 6 cm, and 7 cm. These correspond to the center 
of the source chamber, middle of the double layer, just down- 
stream of the double layer and middle of the downstream 
region. 



proximately 1 x 10 4 m/s. In the next section it will be 
shown that T e ps 4 eV downstream, so this beam travels 
at s» 3c s . This agrees with the expected flow speed if the 
double layer potential drop is the floating potential of 
Eq. V t ~ ^elAfohl/Mi = 3.1c s (for argon). One- 
dimensional cuts of the beam distribution are shown in 
Fig. [TTJ The largest ion-neutral cross section at the ion 
beam energies is charge exchange. This expectation is 
corroborated by the data of Figs. [To] and [TH which show 
that ions lost from the beam show up directly as low- 
energy thermal particles. If the collisions were elastic, 
the beam would slow gradually, which does not happen. 
The ion beams shown in Figs. [10] and [TTJ agree with the 
previous simulations and the ~ 3 c s spe ed downstream 
agrees with previous measurements!^ 1 ^ 
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FIG. 12. The i-directed EVDF for a 0.1 mTorr neutral pres- 
sure simulation at (a) the source sheath edge (x = 1.5 cm), 
(b) just upstream of the double layer (x = 4 cm), (c) just 
downstream of the double layer (x = 6 cm), and (d) the ex- 
pansion chamber sheath edge (x = 9 cm). 




-90 -60 -30 30 60 90 
t. [eV] 



-90 -60 -30 30 60 90 
g. |eV| 



(c) 


* x = 6 cm 







(d) 



x = y cm 




-e"4 
lel2 
1e10 
M 



-90 -60 -30 30 B0 90 -90 -60 -30 30 60 90 

[«V| £. |eV| 



FIG. 13. The x-directed EVDF for a 1 mTorr neutral pressure 
simulation at (a) the source sheath edge (x = 1.5 cm), (b) just 
upstream of the double layer (x — 4 cm), (c) just downstream 
of the double layer (x = 6 cm), and (d) the expansion chamber 
sheath edge {x = 9 cm). 



B. EVDFs and electron temperature 



The EVDF in the x direction is shown in Figs. 12 and 
[T3] for neutral pressures of 0.1 and 1 mTorr. In each 
figure, the EVDF is shown at four positions in the sim- 
ulation domain: the source region sheath edge (s2=1.5 
cm), just upstream of the double layer (DL2=4 cm), just 
downstream of the double layer (DL1=6 cm), and the 
expansion region sheath edge (sl=9 cm). These figures 
can be compared with the model predictions from Fig. [6] 
of Sec. HTDl 

For low neutral pressure (0.1 mTorr), each of the fea- 
tures predicted in Fig. [6] of Sec. |II D| can be seen in the 
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simulation data of Fig. [12] Here, the potential drop of 
the source sheath is |A0 s2 | = 23 V, the double layer 
is |A(/> DL | = 18 V, and the expansion region sheath is 
| A0 s i | = 40 V. At the source region sheath edge, the 
distribution is depleted from the nominal Maxwellian for 
£ x > e|A(/) s | by more than two orders of magnitude. This 
is the truncation due to electron loss to the source bound- 
ary that was predicted in Sec. |IID| The EVDF is also de- 
pleted for £ x < — e|A0 s2 | due to the same wall losses, but 
it has been partially repleted due to scattering over the 
whole the simulation domain. Figure [12] also shows ad- 
ditional depletion for £ x < — e(|A^ s i| + |A</>dl|) = 58 eV 
due to losses to the expansion chamber boundary. Like- 
wise, the predicted features of the EVDF at each of the 
other positions (x = DL2,DL1, and si) compare well 
with the predictions from Fig. [6] 

As the neutral pressure is increased, Fig. |T3] shows that 
repletion of the velocity space intervals subject to wall 
losses also increases. This is simply due to the increase in 
electron-neutral scattering that occurs for higher neutral 
density. The dominant scattering processes for electrons 
on the tail of the Maxwellian (beyond the sheath energy) 
are elastic and ionization collisions. The elastic processes 
cause incident electrons to change velocity by a small 
amount during each scattering event. Repletion of the 
tail happens from a combination of high energy electrons 
scattering from the perpendicular to parallel direction 
and electrons in the parallel direction gaining energy from 
several scattering events. 

Figure [14] shows the +x direction of the EVDF from 
Fig. [13] at positions DL2 and DLL Three populations 
of electrons are present. These include the trapped elec- 
trons below the break energy and tail electrons past the 
break energy that were included in the models of Sec. [TTJ 
The step from one population to the other, which was 
assumed to be a sharp step in the model, is broadened 
due to scattering. Electrons in this intermediate velocity- 
space interval form a third population. Figure [14] also 
shows the effective temperature of each of these three 
intervals. These effective temperatures are calculated 
using a linear least squares fit to the data in the form 
ln(/// Q ) = A£ x + B, in which / is the simulation data 
for the EVDF in the x direction, f a = f(v x = 0) and 
A and B are the constants determined from the lin- 
ear least squares fit. Assuming each interval is close to 
Maxwellian, i.e., straight lines in Fig. [14] the effective 
temperature for that interval is T = 1/\A\. Although the 
step between the trapped and tail populations is not im- 
mediate, the temperature characterizing this interval is 
much colder than either the trapped or tail temperature. 
The model of Sec. [TT| effectively assumes T- mt = 0. 

In the models of Sec. [TTJ it was assumed that both the 
trapped and tail populations had the same effective tem- 
perature T e . Figure [14] suggests that this is a reasonable 
assumption. However, the intermediate population was 
not included in the model and presents a complication in 
that these electrons are effectively colder. In particular, 
we want to determine what temperature should be used 
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FIG. 14. EVDF in the +x direction at the upstream edge 
(x = DL2) and downstream edge (x = DL1) of the double 
layer. Shown are linear least squares fits to three intervals of 
velocity space: trapped, tail and intermediate. Also shown are 
the effective temperatures of each interval. This simulation 
was run with a 1 mTorr neutral pressure. 



in calculating the double layer potential drop. Upstream, 
most electrons are trapped so we expect the total tem- 
perature there to be approximately the temperature of 
the trapped population. However, most of the trapped 
electrons do not contribute to the current balance (see 
Sec. II B). Only electrons that have enough energy to 
escape the double layer £ x > |A</>dl|, i.e., those that 
make it downstream, contribute. Thus, we expect that 



the appropriate temperature to use in calculating the 
double layer potential drop should be the downstream 
temperature. This includes a small part of the trapped 
population [0 < £ x < e(| A<^> s2 | — |A0dl|)]j the whole 
intermediate population, and the whole tail population 
[£ x > e(| A0 S 2| — |A</>dl|)]- As long as the neutral pres- 
sure is within the range that a double layer can form, 
the double layer potential drop is expected to be approx- 
imately the floating potential in which the temperature 
is the downstream temperature: 



|A0 



DL 



2c 



1 + ln 



Mi 
2irm. 



(20) 



Figure [15] shows the electron temperature through- 
out the simulation domain. This is calculated from the 
EVDF using the moment definition: 



T^\ m /^(v-V) 2 / 
3 n J 



(21) 



in which V = i J d 3 v v / is the fluid flow velocity. Also 
shown are three characteristic temperatures of the EVDF 
in each Cartesian direction. Along x, this directional 
temperature is defined as 



Tr 



d 3 v (v x - V x f f 



(22) 



with analogous definitions for the y and z directions. The 
total temperature can be expressed in terms of the direc- 
tional temperatures with the relation: T = (T x + T y + 
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FIG. 15. Electron temperature (T e ) calculated from the 
EVDF using the moment definition of Eq. ( |21[ | (black solid 
line). Also shown are effective temperatures in the x direction 
(green dashed line), y direction (red dash-dotted line) and z 
direction (blue dotted line) calculated using Eq. \22\. 



T 2 )/3. Figure 15 shows that electrons are significantly 



colder in the downstream region than the upstream re- 
gion. This is because the colder intermediate population 
is a greater fraction of the total electron density down- 
stream than upstream. Upstream, most of the electrons 
are in the trapped interval (— e|A0DL| < £x < e|A0DL|) 
and these electrons set the upstream temperature; see 
Fig. [15] Electrons in the x direction are colder than ei- 
ther of the perpendicular directions because the predom- 
inant sink for electron energy is wall losses, which only 
happens in the x direction. Electrons are hottest in the y 
direction because this is the only direction that electrons 
are heated. Figure [15] also shows that there is some elec- 
tron heating from the presheaths of the upstream sheath 
and double layer. Using the 4 eV downstream electron 
temperature from Fig. [15] Eq. pp] ) predicts |A0 DL | w 21 
eV. This agrees well with the approximately 20 eV po- 
tential drop shown in Fig. [8] 



C. Neutral pressure limits 

The potential profile through the simulation domain is 
shown in Fig. [TBJfor neutral pressures of 0.06, 0.1, 2 and 
6 mTorr. The potential drops |A0 s2 |, |A0dl| and | A0 s i | 
are also shown in Fig. [IT] for several neutral pressures 
ranging from 0.04 to 10 mTorr. These were calculated us- 
ing |A0 s2 | = 02 - <j>sw, |A0 DL | = 2 - 01, and |A0 s i| = 
0i — O where <f> sw = <p(x = 0), 02 = <fi(x = 2.5cm), 
0i = 4>{x = 7.5 cm) and O = 0(x = 10 cm) = 0. The 
figures show that as the neutral pressure is decreased, the 
downstream sheath drop increases. The upstream sheath 
and double layer potential remain nearly constant. Sim- 
ulations were also run at 0.01 and 0.02 mTorr, but no 
double layer was found. In these cases the plasma density 
was very low, even though it was stable in time, which 
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FIG. 16. Electrostatic potential profile through the simula- 
tion domain for various neutral pressures. 



is characteristic of there not being enough ionization to 
sustain the discharge. Thus, the minimum pressure to 
sustain the discharge in the simulation was in the range 
between 0.02 and 0.04 mTorr. Although, Fig. [17] shows 
that at 0.04 mTorr the downstream sheath potential drop 
becomes very large (the data point is at 114 V, which is 
off of the figure) and this does not seem physically reason- 
able. Thus, maybe the minimum neutral pressure should 
be considered close to 0.04 mTorr. 



Figure 16 shows that 
for a neutral pressure of 6 mTorr, the potential profile in 
the downstream region is no longer flat, but linearly de- 
creases from the double layer to the downstream sheath. 
This is characteristic of a non- neutral downstream region, 
and the breakdown of the current-free double layer solu- 
tion. At 2 mTorr, the potential in the downstream region 
is flat between the double layer and sheath, suggesting 
that the double layer solution breaks down between 2 and 
6 mTorr. Data points for |02 — 0i| in this high pressure 
region, which are not considered double layer solutions, 



are shown as stars in Fig. 



17 



Equations (11) and (17) provided predictions for the 
minimum and maximum neutral pressures that can sup- 
port a current-free double layer solution. The source and 
downstream simulation domain lengths are L s = Ld = 5 
cm and fo r thermal (« 4 eV) electrons, cr e _ n ~ 1 x 10~ 19 
m -3 |43i4|] Using these parameters, Eq. ( 11 ) yields p mul ~ 
0.03 mTorr and Eq. (17) yields p max ~ 6 mTorr. Both of 



these estimates are consistent with the simulation results 
of Pmin — 0.04 mTorr and p m ax — 2 — 6 mTorr. 



V. SUMMARY 

A model for the EVDF in an expanding plasma with 
a current-free double layer was developed and shown to 
compare well with results of a PIC simulation. The dom- 
inant mechanisms determining the EVDF are depletion 
of high energy electrons due to boundary losses and re- 
pletion of these energy intervals due to scattering. The 
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10"' 10" 
Pressure [mTorr] 

FIG. 17. Upstream sheath potential drop |A0 8 2| (red trian- 
gles), double layer potential drop |A0dl| (black circles) and 
downstream sheath potential drop [A$> s i| (blue squares) as a 
function of neutral pressure. Stars show points for \<f)2 — 4>i\, 
but where the downstream electric field is sufficiently strong 
that it is not considered a double layer solution. 



degree to which these velocity intervals are repleted was 
shown to depend on the ratio of the electron-neutral col- 
lision length to the system size: A e ~ n /L. Assuming a 
simple linear dependence on this parameter, a model for 
the range of neutral pressures that can support a double 
layer was developed. The pressure minimum [Eq. (Ill] 
is determined by the minimum scattering needed to sus- 
tain the discharge. The pressure maximum [Eq. (IT)] is 
determined by current balance through the double layer. 
When the neutral pressure is high, abundant electron 
scattering in the downstream region generates a large 
flux of electrons that can migrate back to the double 
layer and be accelerated by it into the upstream region. 
If too many electrons do this, which happens at high pres- 
sure, current balance across the double layer cannot be 
maintained. The maximum double layer potential drop 
for this configuration is the floating potential using the 
downstream electron temperature. Electrons traveling 
from the downstream to the upstream region causes a 
slight decrease from the maximum. Although this model 
and simulation used a ID domain, the mechanisms of 
depletion due to wall losses and repletion due to scatter- 
ing are expected to be similar in the experiments. These 
results provide information about the EVDF that is es- 
sential for the development of a comprehensive analytic 
model of the experiments. 
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